今天又跳tone啦，想說再把向日葵八號衛星觀測資料及怎麼視覺化說一下。
可以由氣象局開放資料網頁的觀測資料，看到向日葵八號衛星觀測資料有提供很多「頻道」(Band)。

但向日葵八號衛星原始資料有更多頻道，這邊可以看到相關資訊。
如網頁的資訊，向日葵八號衛星總共可提供16個頻道資料，而中央氣象局僅提供其中幾個。

那頻道到底是什麼？只是方便定義而已，真正重要的是頻道對應的波長，我們可以把衛星想像成照相機，地球想像成物件，由於不同的物件對於不同的波長會有吸收或反射等狀況(跟物件自身的電磁波有關)，那相機可以捕捉這些狀態。可以想想看，在黑夜中，如果不開閃光燈，相機是不是什麼都拍不到？其實不是什麼都拍不到，而是所有的光都被吸收，沒有任何光反射至相機的鏡頭給捕捉。這邊有科普介紹。

向日葵八號衛星總共有16個頻道，那一定有其特別的用意，舉個例而言，頻道1至3分別對應的波長剛好是可見光藍、綠及紅的可見光波長。所以如果把頻道1至3透過一些演算法結合後，就可以視覺化可見光的地球，更厲害的話，就可以向氣象局官網提供的真實色地球。

那，氣象局提供的這幾個頻道能夠做什麼呢？這邊有一篇文章主要介紹哪些頻道整合起來後主要的用途。譬如把頻道13與頻道15之差值當成可見光紅色波段、頻道7與頻道13之差值當成可見光綠色波段、頻道13的inverse(就是黑變白，白變黑)當成可見光藍色波段，將之結合視覺化的用途是可以凸顯夜晚微物理物體(雲滴、冰晶等)的資訊。剛好氣象局開放資料有提供可以完成這項產品的所有頻道資料，讓我們試試吧！

我把會用的資料都放在python執行檔該層目錄的day13資料夾內。因為要下載3份資料並讀取，所以就寫了function。

import  requests
import numpy as np
import xmltodict as xdict
from pyhdf.SD import SD, SDC
import cartopy.crs as crs
import matplotlib.pyplot as plt
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature


def downhdf(subdir, fname):
    with open(subdir + "/"+ fname, "r", encoding='utf-8') as fn:
        satdict = xdict.parse(fn.read())
    url = satdict["cwbopendata"]["dataset"]["resource"]["uri"]
    fnout = url.split("/")[-1]
    r = requests.get(url, stream=True)
    with open(subdir + "/" + fnout, 'wb') as fd:
        for chunk in r.iter_content(chunk_size=128):
            fd.write(chunk)

def getright(varhdf):
    scale_factor = varhdf.attributes()["scale_factor"]
    add_offset = varhdf.attributes()["add_offset"]
    hdf_true = varhdf.get()*scale_factor + add_offset
    return hdf_true

for fns in ["O-A0043-002.xml","O-A0043-004.xml","O-A0043-005.xml"]:
    downhdf("day13", fns)

b7 = getright(SD("day13"+ "/O-A0043-002.hdf", SDC.READ).select("temp_3_75um_nom"))
b13 = getright(SD("day13"+ "/O-A0043-004.hdf", SDC.READ).select("temp_10_4um_nom"))
b15 = getright(SD("day13"+ "/O-A0043-005.hdf", SDC.READ).select("temp_12_0um_nom"))

r = np.clip((b13- b15)/255, 0, 1)
g = np.clip((b7 - b13)/255, 0, 1)
b = 1 - np.clip(b13/255, 0, 1)   #inverse

#這邊做一些影響上的處理
gama_r, gama_g, gama_b = 1, 0.5, 0.95
r = np.power(r, 1/gama_r)
g = np.power(g, 1/gama_g) 
b = np.power(b, 1/gama_b)

#把grb結合
rgb = np.dstack([r, g, b])

#畫圖囉，底圖用直角坐標投影
fig = plt.figure(figsize=(12,8))
axs = plt.axes(projection=crs.PlateCarree())
tw_shp = ShapelyFeature(shpreader.Reader("D:\ithome\/tw_shp\COUNTY_MOI_1090820.shp").geometries(), crs.PlateCarree(),facecolor="none",edgecolor='w')
axs.gridlines(draw_labels=True)
ims  = axs.imshow(rgb,extent=(100, 140,9,40),transform=crs.PlateCarree())
axs.add_feature(tw_shp)
filename = SD("day13"+ "/O-A0043-002.hdf", SDC.READ).attributes()["FILENAME"].split("_")
obstime = filename[2] + " " +filename[3] + " UTC"
plt.title(obstime)

畫出來的圖形很尷尬，因為這是看夜晚的，但現在下載資料是在白天阿~~~

晚上再來畫一張試試看

另外補充一下，為什麼我設定imshow的extent範圍是(100,140,9,40)，而不是依照xml檔案提到的經度範圍(92.065 - 140.24)及緯度範圍(8.578 - 45.503)設定為(92.065,140.24,8.578,45.503)。第一個很大的原因是依照xml提供的資訊，畫出來的圖都不正確。所以我就把hdf原始資料的相關資訊print出來，指令如下

print(SD("day13"+ "/O-A0043-002.hdf", SDC.READ).attributes())

上述指令會回傳字典物件，其中有四個key就是我會用(100,140,9,40)的原因，陳列如下
'LAT_SOUTH_SUBSET': 9.0,
'LAT_NORTH_SUBSET': 40.0,
'LON_WEST_SUBSET': 100.0,
'LON_EAST_SUBSET': 140.0,

這就是原因囉。